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Intzoduciion.  We  have  been  working  on  the  following  issues: 

1.  Computations  for  2- Dimensional  and  S-Dimensioual  Inverse  Scattering  Problems 
(ISP)  for  the  Helmholtz  equation.  Convergence  analysis  of  our  numerical  method  for 
this  equation. 

2.  Globally  convergent  numerical  methods  for  3-D  ISP. 

3.  Numerical  method  for  phaseless  1-D  ISP. 

4.  Mathematical  model  for  light  propagation  in  highly  scattering  media 
(interdisciplinary  research  effort). 

Below  we  briefly  describe  these  developments.  They  have  been  described  in  detail 
in  references  [4]- [6],  [8]- [17]. 


1.  Convergence  Analysis  and  Computational  Experiments  for  2-D  and  3-D  ISPs  for  the 
HelmlKdtz  Equation.  See  [4]- [6],  [8]  and  [9]. 

Let  n  be  a  bounded  domain  in  with  a  dielectric  function  l-f£(x),  where 
e(x)  €  L°°(D)  represents  relative  fluctuations  of  the  dielectric  function,  e(x)  =  0  outside 
of  fl  and  £(x)  is  unknown  inside  of  0.  Without  loss  of  generality  we  assume  Jl  C  W, 
where  W  is  the  ball  of  radius  7r/2  with  the  center  at  the  origin.  Let  9W  be  the 
spherical  boundary  of  W  and  S^  be  the  tmit  sphere  in  R^.  If  i/  €  S^,  then  is 

the  scalar  incident  planar  wave  of  the  frequency  k  propagating  in  the  direction  of  the 
vector  V.  Here  <  • ,  •  >  is  the  inner  product  in  R^.  Let  u(x,  v)  be  the  solution  of  the 
forward  scattering  problem 

Au-»-k^(l-f  £(x))u  =  0  (1.1) 

u(x,  I/)  =  -t-u  (x,  u)  (1.2) 

plus  Sommerfield  radiation  conditions. 

The  ISP  Statement.  Let  k  =  const.  >  0  be  fixed.  Determine  the  function  e(x)  in  (1.1) 
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asstoming  that  the  measurement  function  v?{x,  i/),  i/  6  S^,  x  €  &W  is  given  and  the 
function  u(x,  u)  satisfies 

In  [14]  we  proposed  a  quasi-Newton  method  for  this  problem.  During  this  year  our 
developments  in  this  direction  were  the  following: 

a)  We  derived  a  completely  rigorous  proof  of  convergence  of  om  original  algorithm  [4], 

[6]. 

b)  We  tested  complicated  geometries,  rather  than  just  cylinders  as  it  was  earlier  [9]. 
Our  original  code  required  certain  changes  to  do  this. 

c)  The  major  effort  was  devoted  to  solution  of  a  challenging  3-D,  rather  than 
2-D  ISP.  (see  below) 

In  our  algorithm  we  represent  e(x)  in  the  form  of  a  finite  Fourier  series, 

e(x)=  r  ane'<“’''>,  (1.4) 

N<n 


where  N  is  a  fixed  integer,  N  <  ^  in  3-D  case,  and  N  <  v/5k  in  2-D  case.  Hence  we 
cut  high  frequency  harmonics  of  e.  The  algorithm  consists  in  an  iterative  search  of 
Fourier  coefficients  24^. 

Test  ^1.  [9]  As  an  example  of  a  complex  geometry  consider  e(x)  which  is  non-zero 


inside  of  a  non-s3rmmetric  cross  and  a  cylinder,  see  Fig.  la.  That  is 


e(x) 


0.01,  for  ^(x^T^o!27rpTOc2To!2^  <  O.IStt 
0.01,  for  0.6  <  Xj  <  1  and  0.2  <  X2  <  1.3 
0.01,  for  0.2  <  Xj  <  1.3  and  0.6  <  X2  <  1 
0,  otherwise 


We  have  also  included  twelve  randomly  distributed  spikes  in  order  to  see  how  rmdom 


fluctuations  might  affect  our  solution.  Results  of  our  computations  at  N  =  10  are 


displayed  on  Figs.  Ib-ld. 
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Fig.  lb  represents  the  3-D  view  of  the  computed  c;  the  spikes  are  removed.  Figs.  Ic 
and  Id  display  views  from  the  “top”  or  real  and  computed  e  respectively. 
Computations  were  made  on  NAVY  CRAY-Y/MP.  CPU  time  was  2  minutes. 

In  the  test  ^3  we  present  solution  of  a  3-D  ISP,  which  was  obtained,  for  the  first 
time  (cf.  [2],  where  computations  of  a  non-symmetric  3-D  ISP  was  stated  as  a 
challenging  problem).  CPU  time  wm  30  minutes,  at  N  =  7,  which  is  much  less  than 
would  be  required  by  other  currently  available  methods. 

Test  #3.  3-D  ISP. 


In  order  to  avoid  solving  of  3-D  forward  problem,  for  data  simulation  we  have  taken  the 
ball. 


£(x) 


0.02,  for  1x1  <  1.26 
0,  otherwise 


Forward  problem  (1.1),  (1.2)  admits  an  explicit  solution  on  this  case.  However,  we 
have  never  assumed,  whatsoever,  any  kind  of  symmetry  when  inverse  problem  was 
being  solved. 

In  order  to  reduce  dramatically  the  CPU  time,  we  have  made  a  crucial 
modification  of  our  original  algorithm.  Namely  instead  of  using  (2N  -i- 1)^  different 
directions  u  we  employed  some  “basic”  vectors  v.  In  order  to  get  data  for  other  needed 
directions  we  applied  a  special  interpolation  procedure  using  the  fact  that  other  v-s  are 
sufficiently  close  with  basic  ones.  Thus  this  procedmre  allowed  us  to  use  300  “basic” 
directions  u  instead  of  (2  •  7  -f- 1)^  =  3,375  as  it  would  be  the  case  in  [4],  [5]. 

Figs.  2a,  2b  represent  cross-section  of  real  e(x)  by  two  different  pl<mes  at  the 
distance  h  from  the  center  Figs.  2c,  2d  display  corresponding  computed  e.  Finally  on 
Fig.  2e  we  present  cross-section  of  tested  (dashed  line)  and  computed  (solid  line)  e(x) 
by  a  straight  line  through  the  center.  Work  on  an  improvement  of  this  result  is 
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currently  in  progress. 


2.  Globally  Convergent  Numerical  Methods  for  Hyperbolic  3-D  ISP.  [8],  [11],  [12]. 

This  is  an  entirely  new  and  promising  direction  of  the  research  of  numerical 
methods  for  3-D  ISPs.  We  hope  that  these  results  will  lead  us  to  effective  working 
munerical  schemes.  There  are  two  major  {tdvantages  of  these  methods  as  compared 
with  the  method  for  Helmholtz  equation:  (i)  These  are  “global”,  rather  than  “local” 
methods.  That  is  the  basis  of  convergence  is  just  an  aurbitrary  fixed  compact  set,  rather 
than  a  “small”  ball  (in  certain  Banach  space). 

(ii)  These  2ire  “single”,  rather  than  many-source  methods.  That  is  we  consider  non- 
overdetermined  problem  with  just  a  single  source  position. 

Single  source  problems  are  of  much  less  computations  complexity  than  many- 
source  problems  since  they  deal  volumes  of  data  which  are  of  the  order  of  a 
magnitude(8)  less  than  in  the  case  of  many  sources.  The  price  for  these  advantages  is 
two  fold:  (a)  theory  of  these  methods  is  much  more  sophisticated;  (b)  single-source 
problems  are  less  informative. 

The  main  tool  in  this  direction  is  a  deep  modification  of  the  method  of  Carleman 
estimates  [13].  As  to  informativeness  issue,  the  idea  is  to  work  with  the  data  obtained 
from  different  sources  by  a  step-by-step  procedure,  rather  than  working  with  many 
sources  simultaneously.  In  this  case  one  would  improve  the  required  image  step-by-step 
and  CPU  time  will  be  added,  rather  than  mutliplied.  Fin2dly,  one  would  deal  with  a 
much  less  number  of  sources,  which,  in  principle,  should  decrease  experimental 
complexity. 

Let  T  =  constant  >  0,  function  a(x)  6  C^(R^)  and  function  u(x,  t)  is  solution  of  the 
Cauchy  problem 
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=  Au  +  a(x)u,  in  x  (0,  T) 


(2.1) 


u|t^0  =  0>  'itlt=0  =  W 

Consider  cylinder  CL  =  {xj  +  <  r^},  with  the  boundary  u  where  rQ  =  constant  >  0.  Let 

R  =  const.,  R > rQ,  and  n  =  {|x|  <  R}\{CL}. 

The  ISP  statement  Determine  function  a{x)  inside  of  the  domain  D  aissuming  that  a(x) 
is  given  inside  of  CL  and  the  following  function  <^(x,  t)  is  given  as  well 

uL  =  (^(x,  t),  (x,  t)  6  u>  X  (0,  T)  (2.2) 

In  [8],  [11]  we  constructed  a  uniformly  strictly  convex  cost  functional  for  this  ISP.  Thus 
global  convergence  of  a  number  of  numerical  methods  wm  guaranteed.  By  the  method 
[8],  [11]  one  determines  a  finite  nximber  of  Fourier  coefficients  of  a  function  associated 
with  the  wave  field  u(x,  t).  A  more  challenging  and  complicated  problem,  however, 
consists  in  determining  of  a  finite  number  of  Fourier  harmonics  of  the  unknown 
coefficient  a(x)  itself.  During  the  Spring,  Summer  and  Fall  of  1993  we  have  been 
actively  working  on  development  of  such  a  numerical  scheme.  Now  this  method  has 
been  finally  derived  [12].  Thus  algorithm  combines  our  idea  of  quasi- reversibility 
method  [14]  with  the  method  of  [8],  [11]. 

3.  Numerical  Method  for  Phaaeleaa  1-D  ISP.  [17]. 

In  the  classical  formulation  of  ISP  for  1-D  Schrbdinger  equation,  one  has  to 
determine  unknown  potential  V(x)  from  the  reflection  coefficient  R-(k)  [3].  In  many 
important  applications,  however,  only  the  amplitude  |R-(k)p  can  be  measured.  Let 

A  =  {V€L«>(R)nLj(R);  V  is  real 
valued,  V(x)  =  0  for  x  <  0} 
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Here  L2(R)  is  the  weighted  space 

Lj(R)  =  {f:  I  (f(x))(l  +  x2)dx<cx)} 

The  time-reduced  Schrddinger  equation 

^«  +  (k2_V(x)V  =  0,  -oo<x<oo  (3.1) 

then  has  solutions 

$j(x,  k)  =e^^*-HR-(k)e~'^,  x<0 

'w  T(k)e*^,  x->  -I-  oo 

«2(x,  k)  =  T(k)e  -  X  <  0  (3.2) 

~e-i^  +  R  +  (k)e»^,  x^  +  oo 

27ie  ISP  Statement.  Determine  function  V  £  A  assuming  that  the  reflection  amplitude 

r(k)  =  |R-(k)|2,  0<k<oo  (3.3) 

is  given. 

Previously,  we  have  proven  uniqueness  theorem  for  the  problem  (3.1)-(3.3). 
Recently,  we  have  developed  and  tested  a  numerical  method  as  well  [17].  This  method 
works  for  the  case  of  a  relatively  small  number  of  complex  zeros  of  the  function  R(k), 
Fig.  3  represents  exact  (solid  line)  and  computed  (dotted  line)  potentials  V(x). 
Computations  were  made  by  the  algorithm  [17].  Currently,  we  aure  also  working  on 
some  modifications  of  the  idea  of  [16]. 

4.  Mathematical  Model  for  Light  Propagaticm  in  Highly  Scattering  Media.  [10]. 

In  the  past  several  years  different  research  groups  have  been  using  laser-based 
techniques  to  locate  translucent  objects  in  highly  scattering  media  such  as  sea  water, 
biological  tissues,  etc.  [19],  [20].  In  the  future,  this  technique  might  improve  or  even 
replace  conventional  X-ray  tomography.  The  core  of  this  kind  of  imaging,  however, 
must  be  solving  of  full-scale  3-D  ISPs,  since  the  only  hope  to  image  internal  structure  of 
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inclusions  with  a  fine  resolution  lies  in  applying  of  sophisticated  inverse  problems 
numerical  methods.  In  particular,  we  hope  that  our  numerical  methods  might  be 
eventually  applied  to  this  challenging  problem. 

The  first  question  which  comes  in  mind,  however,  consists  in  the  right 
mathematical  model  for  this  kind  of  processes.  Because  of  oxir  already  developed 
algorithms,  it  would  be  desirable  to  have  a  hyperbolic  equation,  which  would  govern 
light  propagation  in  highly  scattering  media.  The  light  emitted  from  ultrafast  laser 
pulses  ( —  100  femtosecond  diuration)  splits  into  three  components  in  highly  scattering 
media;  ballistic,  snake  and  diffuse  [20].  The  ballistic  photons  propagate  along  straight 
lines,  just  as  X-rays.  But  intensity  of  the  ballistic  component  is  so  smeill  that  it  cannot 
be  detected.  So-called  snake  photons  represent  early  arriving  photons.  They  propagate 
slightly  off  straight  lines  and  contain  a  good  portion  of  information  about  inclusions. 
Fortunately,  snake  photons  can  be  detected.  Finally,  diffuse  photons  represent  the  later 
portion  of  the  signal.  They  have  many  scattering  events  before  emission  from  the 
media  and  contain  a  very  little  information  about  inclusions.  Overall,  all  photons, 
except  of  ballistic  ones,  experience  random  scattering  inside  of  such  kind  of  media. 

Commonly,  diffusion  equation  has  been  used  for  description  of  these  processes  [Ij. 
The  immediate  shortcoming  of  diffusion  equation,  however,  is  its  prediction  that 
intensity  of  the  light  will  be  non-zero  instantaneously  everywhere.  Besides,  it  was 
shown  in  [21]  that  diffusion  equation  cannot  describe  snake  photons  and  that  it  is  not 
valid  for  the  thin  media.  After  a  long  and  carefril  analysis  we  suggested  to  use  a  hybrid 
wave-diffusion  equation,  which  is  equivalent  with  the  telegraph  equation  [10].  In  fact, 
wave-type  equations  were  never  used  before  in  this  kind  of  processes,  although  the 
telegraph  equation  was  derived  phenomenologically  for  heat  propagation  in  the  gas  in 
[18,  p.  865]. 
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The  telegraph  equation  for  the  intensity  I(x,  t)  has  the  form 

^Iu  +  ilt-AI  +  ^I  =  0  (4.1) 

We  have  also  proposed  special  initial  conditions  which  reflect  the  collimated  nature  of 
the  initial  laser  beam  propagating  along  the  x^-axis.  These  initial  conditions  are 

I|t,„  =  «(x)  (4.2) 

Here  the  c^6{x.)  term  is  included  to  keep  intensity  from  being  negative  at  t-*oo,  and  ^ 
is  a  small  positive  number  (we  usually  take  ^  —  0.03).  We  can  prove  that  this  term 
provides  a  very  small  impact  until  time  becomes  very  large. 

In  (4.1)  c  is  the  speed  of  light  in  the  media,  D  =  c£j./3,  is  diffusion  coefficient,  is 
the  mean  free  path  of  photons,  and  is  the  absorbtion  length.  In  highly 

scattering  media,  such  as  biological  media,  usually  =  1  -r  3  mm  and  >  100  mm. 

We  have  tested  validity  of  the  equation  (4.1)  by  comparison  of  the  solution  of  the 
problem  (4.1)-(4.3)  with  the  experimental  data  for  a  large  range  of  source- detector 
distances  as  well  as  for  two  source-detector  configurations:  (i)  “face-to-face”,  that  is 
for  the  case  of  transmitted  light;  (ii)  perpendicular.  Experimental  data  were  obtained 
from  Prof.  R.R.  Alfano  (City  College  of  CUNY,  New  York).  Our  data  fitting  results 
show  that  while  for  the  thick  media  both  telegraph  and  diffusion  equation  provide  good 
fit  with  experimental  curves,  for  thin  media  (  <  telegraph  equation  constantly 
provides  much  better  fit. 

Figures  4a,  4b,  5a,  5b,  6a  and  6b  represent  comparison  of  the  solutions  of  diffusion 
and  telegraph  equations  (solid  lines)  with  normalized  experimental  data  (dashed  lines). 
Face-to-face  sotirce-detector  configuration  was  chosen  in  this  particular  set  of 
experiments  (we  also  have  tested  three  more  tests,  and  all  of  them  provided  similar 
results).  Diffusion  equation  under  consideration  has  the  form 
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^It-AI  =  0 

llt=0  = 

Experimental  data  were  arranged  in  such  a  way  that  influences  of  the  boimdaries  of  the 
media  were  negligible.  Figs,  labeled  “a”  and  “b”  represent  fitting  by  the  telegraph  and 
diffusion  equations  respectively.  In  this  set  of  experiments  a  mediiim  with  ^  =  0.55mm 
and  ^2,  =  405  mm  was  taken.  Here  ^  parameters  and  obtained  by  the 

fitting  by  diffusion  equation  for  thick  media.  In  using  telegraph  equation  we  discovered 
that  parameters  and  6.^  should  depend  on  the  source-detector  configiiration.  This  is 
due  to  mixing  up  of  wave  and  diffusion  modes.  Another  explanation  of  this 
phenomenon  is  that  in  non-stationary  kinetics  these  parameters  should  dep>end  on 
photons  momentum.  This  situation  was  predicted  in  [18,  pp.  179-180]. 

Figs.  4a,  4b  jxj/^®  =  18.18,  that  is  the  media  is  thick.  On  Figs.  5a,  b  \x.\f^  =  '1.21\ 
on  Figs.  6a,  b  |x|/^  =  5.45,  and  [xj/l^  =  3.63  on  Figs.  7a,  b.  As  we  clearly  see  from 
these  Figs,  telegraph  equation  provides  a  much  better  fit  for  thin  media  with  |x|/^  <  7. 

One  way  of  improvement  of  this  model  lies  in  working  with  convoluted  data  for 
I(x,  t)  rather  than  with  I(x,  t)  itself.  This  is  especially  true  for  thin  media  because  in 
fact  streak  camera  convolutes  the  data,  and  it  has  10  picosecond  resolution  only. 
Another  way,  might  consist  in  working  directly  with  more  general  transport  equation. 
But  new  numerical  methods  for  ISP  need  to  be  derived  in  this  case,  which  might  open  a 
new  direction  of  our  research.  Overall,  we  have  undertaken  a  significant  effort  in 
theoretical  physics  working  on  this  issue.  Concurrently,  we  consider  possibility  of 
applying  of  our  numerical  methods  for  ISP  to  the  problem  of  imaging  of  inclusions 
hidden  in  highly  scattering  media.  Thus  we  should  work  with  the  telegraph  equation, 
rather  than  just  with  the  wave  equation  (3.1).  Careful  analysis  of  our  algorithms  shows 
that  they  can  be  modified  for  this  case. 
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Figure  Captions 

Fig.  la  Test  ^1  3-D  view  on  the  test  coefficient  e(x). 

Pig.  lb  Test  ^1.  3-D  view  on  the  computed  coefficient  £(x). 

Pig.  Ic  Test  View  from  the  “top”  on  the  test  coefficient  c(x). 

Fig.  Id  Test  ^1  View  from  the  “top”  on  the  computed  coefficient  e(x). 

Fig  .  2a  Test  ^2  Cross-section  of  the  test  coefficient  £:(x)  by  the  phase  at  h  =  0  from 
the  center. 

Fig.  2b  Test  ^2  Cross-section  of  the  test  coefficient  e(x)  by  the  plane  of  h  =  0.26  from 
the  center. 

Fig.  2c  Test  ^2  Cross-section  of  computed  £(x)  by  the  plane  at  h  =  0  from  the  center. 

Fig.  2d  Test  ^2  Cross-section  of  computed  £(x)  by  the  plane  at  h  =  0.26  from  the 

center. 

Fig.  2d  Test  ^2  Cross-section  of  tested  (dashed  line)  and  computed  (solid  line)  £(x) 
by  a  straight  line  through  the  center. 

Fig.  3  Pha^less  1-D  ISP.  Exact  (solid  line)  and  computed  (dashed  line)  potential 

Fig.  4a  Telegraph  equation  solution  (solid  line)  and  experimental  data  (dashed  line) 
for  lxl/«J  =  18.18.  V  /  V  / 

Pig.  4b  Diffusion  equation  solution  (solid  line)  and  experimentad  data  (dashed  line) 

for  |x|/€f  =  18.18.  V  y  v  ; 

Fig.  5a  Telegraph  equation  solution  (solid  line)  amd  experimental  data  (dashed  line) 
for  |x|/IJ  =  7.27. 

Fig.  5b  Diffusion  equation  solution  (solid  line)  and  experimental  data  (daished  line) 
for  |xl/«5  =  7.27. 

Fig.  6a  Telegraph  equation  solution  (solid  line)  and  experimental  data  (dashed  line) 
for  |x|/£j  =  5.45.  r  v  / 

Fig.  6b  Difr^ion  equation  solution  (solid  line)  and  experimental  data  (daished  line)  for 
|x|/£j  =  5.45. 

Fig.  7a  Telegraph  equation  solution  (solid  line)  amd  experimentail  data  (daished  line) 
for  \x\lq  =  3.63. 

Fig.  7b  Diffi^ion  equation  solution  (solid  line)  and  experimental  data  (dashed  line)  for 
l*l/^t  “  3.63. 
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